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ABSTRACT 

In  this  paper  we  present  an  overview  of  classical  results 
about  the  variance  reduction  technique  of  control  variates. 
We  emphasize  aspects  of  the  theory  that  are  of  importance  to 
the  practitioner,  as  well  as  presenting  relevant  applications. 

1  INTRODUCTION 

The  method  of  control  variates  is  one  of  the  most  widely 
used  variance  reduction  techniques.  Its  popularity  rests  on 
the  ease  of  implementation,  the  availability  of  controls,  and 
on  the  straight  intuition  of  the  underlying  theory. 

To  keep  the  presentation  simple,  we  frame  the  results 
for  the  case  where  the  parameter  to  be  estimated  is  a 
scalar  in  the  setting  of  terminating  simulations,  although 
the  theory  extends  to  the  multi-response  setting  and  to  the 
steady-state  simulation  context.  The  emphasis  is  on  creating 
valid  confidence  intervals  and  on  understanding  the  variance 
reduction  achieved  by  the  estimator. 

For  the  remainder  of  Section  1  we  outline  the  paper 
and  discuss  the  relevant  literature.  In  the  second  section 
we  present  the  basic  formulation  of  control  variates,  which 
includes  finding  the  optimal  control  coefficient  and  creating 
an  asymptotically  valid  confidence  interval.  Because  the 
optimal  control  coefficient  is  generally  unknown,  we  discuss 
the  loss  of  variance  reduction  caused  by  its  estimation  and 
introduce  the  idea  of  loss  factor. 

Section  3  presents  the  relationship  between  control 
variates  and  the  method  of  regressions.  This  relationship  is 
useful  to  obtain  an  expression  about  the  limiting  variance 
of  the  control  variate  estimator  that  uses  an  estimate  of 
the  optimal  control  variates  coefficient.  We  also  show  the 
relationship  between  least  squares  regression  and  control 
variates. 

In  Section  4,  we  discuss  the  method  of  batch  means  as  a 
way  to  overcome  the  bias  introduced  in  the  control  variates 
estimator  when  the  optimal  control  variates  coefficient  needs 
to  be  estimated.  We  close  that  section  by  presenting  an 
asymptotically  valid  confidence  interval. 


Section  5  deals  with  non-linear  control  variates.  We 
show  that  these  type  of  control  schemes  are  no  more  efficient, 
in  terms  of  variance  reduction,  than  linear  control  variates. 

In  the  last  section  we  present  some  applications  of 
control  variates  in  the  realm  of  finance.  We  make  use 
of  the  examples  to  illustrate  the  more  general  problem  of 
finding  and  selecting  control  variates. 

The  literature  in  the  theory  and  applications  of  control 
variates  is  quite  extensive,  and  we  do  not  intend  to  provide 
an  exhaustive  list  here.  The  paper  by  Nelson  (1990)  and  the 
work  of  Loh  (1995)  contain  a  very  complete  list  of  relevant 
references.  We  also  recommend  the  paper  of  Lavenberg 
and  Welch  (1981),  which  was  the  first  to  give  a  complete 
and  rigorous  exposition  of  control  variates.  At  a  more 
introductory  level,  the  books  by  Bratley  et  al.  (1987)  and 
Law  and  Kelton  (1991)  provide  the  fundamentals  of  control 
variates. 

2  PROBLEM  FORMULATION 

AND  BASIC  RESULTS 

We  start  by  considering  the  problem  of  estimating  by  Monte 
Carlo  simulation  a  scalar  parameter  a  that  can  be  expressed 
as  the  expectation  of  a  random  variable  Y ,  that  is  a  —  EY. 
Let  Yi  be  an  output  of  the  ;’th  iteration  of  the  simulation, 
done  in  a  way  so  that  the  replications  Y\ ,  Yi. ... .  Yn  obtained 
after  n  iterations  are  independent  and  identically  distributed 
(i.i.d.)  as  the  random  variable  Y.  The  natural  point  estimator 
of  a  is  the  average. 


The  method  of  control  variates  arises  when  the  simu- 
lationist  has  available  a  random  vector  C  e  W1  (notation 
point:  vectors  are  columns  and  ■'  denotes  transpose)  with 
known  mean  /xc  that  is  jointly  distributed  with  Y  as  an 

additional  simulation  output.  Let  Ci,C2 . C„  be  the 

sequence  of  i.i.d.  observed  outputs.  The  idea  is  to  use 
the  “error"  l/n^”=1C;  —  //c  to  control  Yn.  Intuitively, 
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when  Y  and  C  are  positively  correlated,  we  should  intro¬ 
duce  1  / n  Q  —  Me  in  a  manner  so  that  we  adjust  Y„ 

upwards  when  C„  <  /xc  and  downwards  when  C„  >  //  c . 
One  manner  to  achieve  that  is  via  the  linear  transformation, 

Yn(X)  =  Yn  -  XT  ^  C'  -  ^  -  (!) 

where  the  vector  X  e  Rrf  is  chosen  to  minimize  Var  F„(A). 
Assume  that  E{Y2  +  C1  C)  <  oo,  and  let  <ryc  be  the  d- 
dimensional  vector  whose  elements  are  the  covariances  of 
Y  with  each  of  the  d  components  of  C.  Then,  X  is  chosen 
so  that, 

A.  =  argmin  JvarF  —  2Xl  crYC  +  X1  ECC1  A  j  . 

If  we  assume  that  ECC1  is  non-singular,  then  the  first 
and  second  order  optimality  conditions  of  the  minimization 
problem  imply  that  there  exists  a  unique  optimal  solution, 


Student’s-t  with  n  —  1  degrees  of  freedom,  and  we  can 
construct  the  confidence  interval, 

P(Yn(X)  -tn- i(l  -  y/2)  Var(F„(A))1/2  <  a 

<  F„(A)  +tn- i(l  -  y/2)  Var {Yn{X))x'2)  «  1  -  y, 

where  f„_  i  ( 1  —  y/2)  is  the  1  —  y/2  quantile  of  the  t- 
distribution  with  n  —  1  degrees  of  freedom  for  0  <  y  <  1 . 
When  the  random  vector  (F,  C)  is  multivariate  normal  the 
last  approximation  is  exact. 

In  general,  however,  the  covariance  structure  of  the 
random  vector  (Y,  C)  is  not  fully  known  prior  to  the  simu¬ 
lation.  An  efficient  approach  to  overcome  this  difficulty  is  to 
estimate  these  moments  from  the  already  available  samples 
(Yj,  C ;),  i  —  1, . . ,  ,  n.  The  relevant  unbiased  statistics  for 
the  second  moments  are 

1  " 

Sc(n)  = - -  V(C /  -  C„)(C i  -  C„)T, 

n  —  1  z ' 
i=  1 


X  =  (E  CCT)~l(TyC. 


With  this  choice  of  X  the  variance  reduction  achieved 


is 


Var(T„  (A)) 
Var(F„) 


=  1  -  R 


¥  C’ 


(2) 


where  R^c  =  aJc(ECCT)~1crYC/'VarY  is  the  square  of  the 
multiple  correlation  coefficients  between  Y  and  C. 

The  following  assumption  will  hold  throughout  the 
paper. 

Assumption  1.  (Functional  Central  Limit  Theorem) 
Assume  that  the  stochastic  processes  Y  —  (Y (t)  e  R  : 
t  >  0)  and  C  =  (C(f)  e  :  t  >  0)  are  outputs  of  the 
simulation,  and  define  X  =  (X(f)  =  (Y(t),  C(t))  :  t  >  0). 
Let, 

1  fnt 

X„(0  =  -  /  X(s)ds. 
n  J  o 


We  assume  that  there  exists  a  constant  fix  e  Rrf+1  and 
a  positive-definite  matrix  X  e  )  such  that  the 

following  limit  holds: 


n1/2(Xn(t)  -  iixt)  X1/2B(f),  for  0  <  t  <  1, 


as  n  — >•  oo,  where  B(-)  is  a  standard  Brownian  motion  in 
Rd+1,  and  the  convergence  is  weak  in  the  space  Z)[0,  1] 
(a  good  reference  on  this  topic  is  Billingsley  1999).  In 
this  paper  we  study  the  process  Y  ( t )  =  Y\_t j  along  with  the 
vector  valued  control  process  C  it)  =  C  - 

One  of  the  key  issues  of  every  simulation  is  to  assess  the 
accuracy  of  the  final  estimator  via  confidence  intervals.  The 
distribution  of  ( Yn  (A)  —  a)/  Var (Yn  (A))  is  approximately 


and. 


Sy(n)  = 


£(r«  -  Yn)2, 

i=l 


i  » 

Syc(n)  =  - -  Yjyi  -  Yn)(Ci  -  C„). 

n  —  1  ' 

In  terms  of  these  statistics,  the  optimal  control  parameter 
A„  is 

An  =  Sc(n)-]  Syc(n). 

The  modified  point  estimator  for  a  now  becomes 


YniXn)  =  T„  -Aj(C„  -  Mel¬ 


on  e  problem  with  this  approach  is  that  we  cannot 
use  the  t-statistic  to  generate  confidence  intervals  because 
the  controlled  output  replicates  [Ij-  —  Xjt  (C;  —  //t)],  i  — 
1 , ,n,  are,  in  general,  dependent  of  each  other.  Under 
the  additional  assumption  that  (F,  C)  has  a  multivariate 
normal  distribution  in  Rd+1,  however,  it  can  be  shown  (see 
Lavenberg  and  Welch  1981  and  Lavenberg  et  al.  1982) 
that  F„(A„)  is  unbiased  and  an  expression  for  Var(F„(A;!)) 
can  be  obtained  by  standard  regression  techniques.  It  also 
can  be  shown  that  (F„(Ah)  —  a)/  Var(Fn(A„))1/2  has  a  t- 
distribution  with  n—d—l  degrees  of  freedom,  which  permits 
the  creation  of  an  exact  (1  —  y )  1 00%  confidence  interval. 


P(Yn(Xn)  -  tn-d- i(l  -  y/2)  Var(F„(A„))1/2  <  a 
<  F„(A„)  +  tn-d- 1 (1  -  y/2)  Var(F„(A„))1/2)  =  1  -  y. 
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The  variance  reduction  reduction  achieved  in  this  setting 
is  (more  details  in  Lavenberg  and  Welch  1981  and  Lavenberg 
et  al.  1982) 

Var(7„(A„))  _  n-2  2 

- = -  —  - (1  JXyr-'). 

Var(7„)  n-d-2  rL 

The  factor  (n  —  2 )/{n  —  d  —  2)  >  1  determines  the 
variance  increase  due  to  the  estimation  of  the  covariance 
structure  of  (T,  C)  (compare  with  equation  (2)),  and  for  this 
reason  it  is  called  the  loss  factor.  We  say  that  a  sequence  of 
random  vectors  ( Xn  ■  n  >  1)  is  op(n~ l/2)  if  nl/2Xn  =>■  0  as 
n  — >•  oo,  and  that  (/„  :  n  >  1)  is  o(a„)  a.s.  if  Xn/an  — >■  0 
a.s.  as  n  —>  oo. 

The  following  expansion  of  the  loss  factor, 

Var(7„(A„))  k 

- = - =  1  H - Ion  ), 

Var(T„  (A,))  n 

asserts  that  the  variance  reduction  loss  caused  by  estimating 
A  converges  to  zero  at  a  rate  d /n.  In  addition,  for  any  con¬ 
sistent  estimator  A„  of  A,  under  further  uniform  integrability 
assumptions  (see  Loh  1995,  p.  22)  and  Assumption  1,  the 
following  result  holds: 


Var(7„(A„)) 

lim  - = - =  1. 

«^oo  Var(T„  (A)) 


(3) 


Equation  (3)  asserts  that  as  long  as  the  estimator  of  A 
is  consistent,  we  will  obtain  the  1  —  Ryc  variance  reduction 
guaranteed  by  7„(A). 

The  results  of  this  section  can  be  extended  to  the 
multiresponse  setting,  in  which  case  7  and  a  are  vectors  in 
R9,  q  >  1;  see  Rubinstein  and  Marcus  (1985)  for  details. 


3  THE  REGRESSION  APPROACH 
TO  CONTROL  VARIATES 


We  mentioned  in  the  last  section  that  an  expression  for 
Var(7„(A„))  can  be  obtained  using  regression  techniques. 
In  this  section  we  make  the  connection  between  regression 
analysis  and  control  variates  in  the  bivariate  normal  setting 
(with  just  one  control  variable)  more  explicit. 

Specifically,  assume  that  the  pairs 
(Ti,  Ci),  (Y2,  C2), . . . ,  (T„,  C„)  are  i.i.d.  bi¬ 

variate  normal.  Conditional  on  C;  we  have 
E(Yi\Ci)  =  a+A(C;  —  pc)  where  A  =  —Co  v(7,  C)/VarC, 
and  Var(7,jC,)  =  VarT(l  —  pyc),  where  pyc  is  the 
correlation  coefficient  between  Y  and  C.  Then  we  can 
express. 


where  the  e-s  are  i.i.d.  normal  with  mean  zero  and  vari¬ 
ance  Var  7 ( 1  —  PyC),  and  independent  of  (7,-,C()  for 
i  —  1 , ,n. 

The  minimum  least  squares  problem  is  to  find  param¬ 
eters  a  and  A  that 


minimizeQ.  ^  ef. 

i=  1 

Solving  the  problem  above  we  obtain  a  —  Y„  —  A (C„  —  pc) 
and  A  =  Syc(n)/Sc(n),  the  least  squares  estimators  of  a 
and  A  respectively.  The  variance  of  a  can  be  found  to  equal 
(see  Lavenberg  and  Welch  1981) 

Vara  =  S2(n){n~l  +  (n  -  1)-1(C„  -  p c)2/Sc(n )), 

where  S2(n)  —  ^f^(Sy(n)  —  S2c(n)/Sc(n)).  This  result 
can  be  extended  to  the  multi-control  case  to  provide  an 
expression  for  Var(  Yn  (A„ )). 

Relaxing  the  normality  assumption,  getting  the  best 
linear  fit  of  the  pairs  (Fj,  Ci),  (Y2,  C2), . . . ,  ( Yn ,  C„ )  by 
minimizing  ^"=1  (Y,  —  a  —  k(Cj  —  iic))2  in  a  and  A  yields 
a  =  Yn  —  A (C„  —  iic)  and  A  =  Syc(n)/Sc(n )  respectively. 
That  is,  the  line  a  —  A (C„  —  pLc)  passes  through  the  points 
( Yn ,  C„)  and  (a,  /xf),  so  that  Yn  is  adjusted  to  a. 

4  BATCH  MEANS 

As  already  mentioned  in  Section  1,  a  problem  with  us¬ 
ing  the  estimated  optimal  coefficient  A„  is  that  the  con¬ 
trolled  simulation  outputs  (K|  —  A„  (C 1  —  pc),  Y2  ~  A „  (C2  — 
Pc),  . ...  Yn  —  A„ (C„  —  He))  are  no  longer  independent  of 
each  other.  This  precludes  the  computation  of  confidence 
intervals  for  a,  and  Y„  (A„)  is  no  longer  unbiased.  The  same 
issue  happens  in  the  steady-state  simulation  context.  The 
batch  means  method  tackles  this  problem  by  splitting  the 
output  into  a  fixed  number  of  m  batches  (with  mk  =  n, 
the  integer  k  being  the  batch  size),  and  forming  a  batch 
means  sequence  with  elements  that  are  asymptotically  i.i.d. 
normal  as  n  — »  00.  As  in  Section  2,  our  main  objective 
is  to  determine  a  valid  confidence  interval,  and  to  find  the 
variance  reduction  achieved  by  this  method. 

More  precisely,  let  X,  =  (7,-,  C;),  for  i  —  1 
For  the  batches  i  —  1,  . . . ,  m,  define  the  batch  means  by, 

1  ik 

*(i,n)=-  J2  XJ> 

j=d- 1)*+1 


Yj  —  a  +  A  (Ci  —  pc)  +  €{,  for  i  =  1, . . . ,  n, 


with  7 (i,  n)  (resp.,  C (i,  n ))  being  the  first  (resp.,  second  to 
(d  +  l)’th)  component(s)  of  X(i,  n ).  We  argue  by  induction 
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that  the  elements  (X(i,  n)  :  i  —  1 ,  m)  are  asymptotically 
i.i.d.  normal.  Consider  the  first  batch, 

X(l,  n)  =  mX„(l/m). 

By  Assumption  1, 

/  1  A1/2 

mX„(l/m)  =>•  (  —  £j  N(nx,I), 

where  £  is  the  matrix  with  first  row  (Var  T,  n  jc)  and  second 
row  (ervc,  ECCt),  and  I  is  the  d  x  d  identity  matrix. 
Considering  now  X((7  +  1),  n)  we  have, 

X((f  +  1),  n)  =  m(Xn((i  +  1  )/m)  -  X„(//m)). 

Using  Assumption  1, 


We  want  to  compare  the  variance  reduction  achieved  by  this 
estimator  with  that  of  the  batch  means  controlled  estimator 
in  terms  of  X ,  which  is 

1  m  i  m 

iv„a)  =  -  V  f(i,  n)  -  A/-  V(C(«,  n)  -  Me), 
m  z — '  m  z — ' 

i=l  i=l 

Since  the  batch  means  X(i,  n)  are  asymptotically  i.i.d. 
normal  as  n  ->  oo,  the  number  of  batches  m  in  the  equa¬ 
tion  above  is  the  analogous  of  the  number  of  replications 
n  of  Equation  (1).  Indeed,  under  certain  uniform  inte- 
grability  conditions  on  the  sequences  (Ymn(X(m,  n)))„  and 
(Ym  n (X))n,  and  Assumption  1,  one  can  show  (see  Loh  1995, 
p.  37)  that  for  m  >  d  +  2, 

Var  Ym  n(X(m,  n))  m  —  2 

Var  Y„un(X)  m—d  —  2' 


m (Xn((i  +  1  )/m)  -  X„  (i/m))  -  Mx 


.1/2 


B 


(\  \ 1/2 

(-Sj  N (0,  /). 


The  Brownian  motion  term  above  is,  by  the  independent 
increments  property  of  Brownian  motion,  independent  of 
v1/2w  / m )  —  B((i  —  1  )/m))  (which  is,  by  the  induction 
hypothesis,  the  asymptotic  distribution  of  X(i,  n)).  The 
asymptotic  normality  of  the  batch  means  ensures  that  we 
can  construct  confidence  intervals  that  are  asymptotically 
valid. 

The  relevant  statistics  in  the  batch  means  context  are 


as  n  — >  oo.  Consequently  there  is  a  tangible  loss  of  variance 
reduction  for  the  batch  means  estimator  that  uses  an  estimate 
of  the  optimal  control  coefficient  X. 

One  of  the  advantages  of  the  batch  means  approach  is 
that  it  allows  the  creation  of  asymptotically  exact  confidence 
intervals.  It  can  be  argued  (see  Nelson  1990  and  Loh  1995) 
that, 

Y,nj,(X{m ,  »))  -  a 
Var 

converges  in  distribution  to  a  Student’s-t  random  variable 
with  m  —  d  —  1  degrees  of  freedom  as  n  —*■  oo.  So  the 
batch  means  method  decreases  the  number  of  degrees  of 
freedom.  A  confidence  interval  can  be  generated. 


1  ' 

Sc(m ,  n )  = - -  V/C (i,  n)  -  C„)(C (i,  n)  -  C„)T. 

m  —  1  L — 4 
i= 1 

and, 

^  m 

Sycim,  n)  = - -  T(Y (i,  n)  -  y„)(C(i,  n)  -  C„). 

m  —  1  L — 4 
i=  1 

Letting  X(m,n)  —  Sc(m,  n)_1  Syc(m,  n),  the  batch 
means  controlled  estimator  in  terms  of  X(m ,  n)  is 


P(Ym,n(X(m,  n)) 

-tm-d- 1 (1  -  y/2)  Var(fm,„(A.(m,n)))1/2  <a 
<  Ym  rl(X(in,  n )) 

+  fm-<f-t(l  -  y/2)  Var(fm,„(A.(w,  n)))1/2)  ->  1  -  y, 
as  n  — »•  oo. 

Our  presentation  of  the  method  of  batch  means  makes 
clear  that  selecting  the  appropriate  number  of  batches  m 
is  an  important  decision  for  the  analyst;  this  issue  is  well 
explained  in  Nelson  (1990). 


j  m 

Ym,n(X(m,n ))  =  —  K (i,  n) 


1=1 


—  X(m,  n)T —  (C (i,  n)  —  Me)- 


i=l 


5  NON-LINEAR  CONTROL  VARIATES 

In  this  section  we  consider  the  performance  of  control 
variates  when  they  are  related  to  Yn  in  a  non-linear  way. 
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The  results  presented  in  this  section  are  contained  in  Glynn 
and  Whitt  (1989).  For  example, 

-  Cn 
F„  — , 

Me 


would  be  two  such  schemes  when  Cel.  More  generally, 
we  deal  with  a  scalar  function  /  with  domain  in  Rrf+1.  The 
function  /  of  the  last  two  examples  is  f(y,  c )  =  yc/ Me  and 
f(y,c )  =  ycnl<-  respectively,  and  satisfies  f(y,  /xf)  =  y. 
This  last  property  ensures  that  f(Y„,  C„)  =>-  or  if  ( Yn ,  C„)  =>■ 
(a,  Me),  so  we  only  will  consider  such  functions  in  the 
discussion  that  follows. 

The  variance  reduction  associated  with  any  given 
function  /  will  depend  on  the  limiting  variance  of 
n1/2(/(f„,  C„)  —  a).  Now,  when  /  has  continuous  first 
partial  derivatives  in  a  neighborhood  around  (a,  Me),  we 
can  obtain  via  Taylor’s  theorem  a  first-order  linear  approx¬ 
imation  of  /  around  f(a,  /xc), 

f(Yn,  C„)  =  / (a,  Me)  +  (Yn  -  a ,  C„  -  Me )rV/(f„,  e„), 

where  the  random  variable  §„  and  the  random  vector  e„  e  W1 
lie  on  a  segment  with  end-points  ( Yn,a )  and  (C„,Mc) 
respectively.  Since  (F„,C„)  =>■  (a,  Me),  we  also  have 
(f„,  e„)  =>■  (a,  Me),  and  we  can  write 

f(Yn,  C„)  =  a  +  (F„  -  a,  C„  -  Mc)rV/(a.  Me) 

+  Op(n“1/2). 

Note  that  V/(a,  /xc)  =  1,  so  that, 

f(Yn,  C„)  =  f„  +  (C„-Mc)rVc/(a,  /xc)+Op(n_1/2),  (4) 

where  Vc/  is  the  vector  of  partial  derivatives  of  /  with 
respect  to  the  C  components.  Thus,  the  limiting  distribution 
of  nl'2f(Yn,Cn)  is  the  same  as  that  of  the  linear  control 
nl'2(Y„  —  X1  (Cn  —  Me)),  with  Vc/(a,  fic)  standing  in  lieu 
of  —  X.  This  result  implies  that  non-linear  control  variates 
cannot  improve  the  variance  reduction  achieved  by  linear 
control  variates,  in  the  limit  as  n  —>  oo.  Indeed,  equation 
(4)  results  in, 

rc1/2(/(f„,C„)-a) 

=  nx/2(F„-«) 

+  n1/2(C„  -  Mc)rVc/(a,  Me)  +  op(l). 


Sending  n  — >  oo.  Assumption  1  and  the  converging-together 
lemma  imply  that, 

nl'2(f(Yn,Cn)  -  a)  =►  N(0,aj), 

where, 

a2  =  Var  F  +  2 Vc/(a,  Me)7  o-.vc 

+  Vc/(ar,  Mc)r£(CCr)Vc/(«,  Me)- 

Selecting  /  so  that  Vc/(a,  Me)  =  —  (£’CCr)_1cr),c 
is,  according  to  the  discussion  of  Section  2,  the  variance 
minimizing  function.  This  selection,  in  turn,  implies  that 
the  optimal  control  variate  function  is  linear  with  control 
coefficient  given  by  —  (ECC1  )^1cryc. 

6  APPLICATIONS  OF  CONTROL  VARIATES 

In  this  section  we  present  several  applications  of  control 
variates  in  finance.  Our  first  example  uses  what  are  called 
“internal"  control  variables,  so  called  because  the  control 
variables  are  random  variables,  or  functions  of  them,  used 
as  an  input  to  the  simulation  model  and  are  often  easy 
to  parameterize.  One  of  the  advantages  of  using  internal 
control  variables  is  that  the  additional  computational  cost 
incurred  by  adding  them  is  usually  small  relative  to  the 
overall  cost.  Normal  random  variables  are  often  used  in 
finance  to  drive  pricing  models,  suggesting  the  use  of  their 
known  mean  and  variance  as  control  variables. 

Another  example  of  internal  control  variables  (bor¬ 
rowed  from  Szechtman  and  Glynn  2001)  is  provided  by  the 
computation  of  an  Asian  option  via  simulation  under  the 
risk-neutral  measure  (see  Duffle  1996  for  details).  More 
specifically,  assume  that  the  price  process  §  =  (f(s)  :  s  >  0) 
of  the  underlying  asset  is  geometric  Brownian  motion,  and 
let  k  denote  the  strike  price.  Then,  the  price  of  the  Asian 
option  is  given  by  the  expectation  of  the  random  variable 
F  given  by. 


where  (a)+  =  max(a,  0)  for  a  scalar  a. 

In  this  context,  one  can  analytically  find  the  expectation 
of  the  integrand  in  the  last  equation, 

C=  f  i;(s)ds, 

Jo 

to  be  EC  —  2(exp(l/2)  —  1).  Therefore  C  can  be  used  as 
an  internal  control  for  F. 

External  controls  are  controls  that  are  jointly  distributed 
with  the  replicates  of  the  random  variable  whose  expectation 
we  wish  to  estimate,  and  that  are  generated  in  addition  (often 


148 


Szechtman 


as  a  separate  model  driven  by  the  same  random  input  as 
the  main  model)  to  the  main  model. 

From  Glasserman  (2003),  we  show  an  example  of  ex¬ 
ternal  control  variables.  We  want  to  price  an  option  with 
expiration  time  T  with  strike  price  k  and  whose  underlying 
asset  price  5(f)  has  dynamics  driven  by, 

dS(t) 

- =  rdt  +  <r(t)d  B(t), 

5(f) 

where  the  volatility  a(t)  may  be  random  or  a  function  of 
S(t).  In  order  to  simulate  the  price  dynamics,  we  simulate 
5  at  discrete  times  t\, ...  ,tn  —  T  via  the  recursion, 

Siti)  9 

— - -  =  exp([r  -  1/2ct (f;_i)2](f;  -  ti-i) 

j(h'-l) 

+  a(ti-i)(ti  -  ti-i)1/2Zi), 

where  the  Z,  ’s  are  i.i.d.  standard  normal  random  variables 
and  a(ti)  is  driven  by  its  own  recursion.  The  idea  is  to  run 
another  simulation  alongside  with  constant  volatility  a  and 
initial  condition  5(0)  =  5(0), 

S(tj) 

S(ti- 1) 

=  exp([r  -  1/2(7 2](f;  -  ti- 1)  +  o(ti  -  f;_i)1/2Z;), 

where  the  Z,  ’s  are  the  same  (common  random  numbers) 
as  in  the  model  for  5.  If  the  price  of  the  underlying  asset 
follows  a  geometric  Brownian  motion,  then  we  can  use  the 
Black-Scholes  formula  to  find  E(S(tn )  —  k)+  analytically. 
With  these  assumptions,  generate  controlled  replications, 

(5(G)  -k)+-k  ((5(G)  -  k)+  -  £(5(G)  -  k)+)  , 

to  form  the  usual  control  variates  estimator. 

The  efficiency  of  this  approach  will  depend,  among  other 
factors,  on  judiciously  choosing  the  constant  volatility  <7. 
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